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We compare the semiclassical Hartree-Fock approximation for a trapped Bose gas to a direct 
Path Integral Quantum Monte Carlo simulation. The chosen parameters correspond to current *^Rb 
experiments. We observe corrections to the mean-field density profile. The Path Integral calculation 
reveals an increase of the number of condensed particles, which is of the same order as a previously 
computed result for a homogeneous system. We discuss the experimental observability of the effect 
and propose a method to analyze data of in-situ experiments. 
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The experimental realization of Bose-Einstein conden- 
sation (BEC) in dilute atomic vapors lias gener- 
ated extraordinary experimental and theoretical interest. 
From the viewpoint of quantum many-body physics, the 
trapped atomic vapors are peculiar. Well above the crit- 
ical point, the gases are extremely dilute, and their de- 
scription as noninteracting bosons is very accurate. As 
the condensation sets in, the trapped atoms are strongly 
compressed in real space. Then, interactions become 
much more important and even the simplest thermody- 
namic quantities (spatial distribution, condensate frac- 
tion, etc.) have to be obtained by the appropriate quan- 
tum many-body technique. At zero temperature the Bo- 
goliubov approach of weakly interacting Bose gases is well 
established [Q. There, the macroscopic condensate wave 
function is given by the Gross-Pitaevskii equation 
The effects of noncondensed particles at finite tempera- 
tures can be included via the Hartree-Fock-Bogoliubov 
equations, which have been solved in the Popov approx- 
imation ^ and in various simplified forms 0. 

The equilibrium properties of Bose gases can also be di- 
rectly computed by Path Integral Quantum Monte Carlo 
(QMC) simulation without any essential approximation 
1^. The QMC approach was used to first show ||^ that 
the trapped Bose gas with repulsive interactions has a 
lower critical temperature Tc and a smaller number A^o 
of condensed particles below Tc than the free Bose gas, 
one of the predictions of mean- field theories 

Very importantly, the QMC calculation is free of sys- 
tematic errors: it gives a somewhat noisy, but otherwise 
exact numerical solution of BEC. For simulations corre- 
sponding to dilute atomic vapors, the QMC calculation 
can be performed directly for the large particle numbers 
10^ — 10^) and the temperatures of experimental in- 



terest. In this situation, there is little need to perform 
finite-size scaling. QMC provides us thus with a unique 
opportunity to check the above-mentioned many-body 
techniques directly in an experimental setting, a com- 
parison which we provide here. 

It is the main point of this paper to study the correc- 
tions of a full many-particle treatment to a mean-field 
theory. We compare the QMC calculations to the semi- 
classical Hartree-Fock (HF) approximation, which deter- 
mines the local density in accordance with the local trap 
potential and the mean interaction energy between the 
particles. The HF approximation provides an ideal ref- 
erence system, since it corresponds locally to a homoge- 
neous ideal Bose gas. The corrections to HF should be 
therefore comparable to previous calculations in a homo- 
geneous system JT2|-p^. 

For a system closely corresponding to many experi- 
ments with ^^Rb we find deviations from the HF density 
which are concentrated in the overlap region between the 
excited atoms and the condensate. The number of con- 
densed particles Nq is increased around the critical tem- 
perature by about 5%. We also discuss a scheme for a 
detailed analysis of experimental in-situ measurements. 

The QMC calculation is beset with statistical fluctua- 
tions. In fact, the noise of independent Monte Carlo con- 
figurations reproduces the sample-to-sample variations of 
repeated experimental measurements at the same tem- 
perature T and particle number N (or chemical poten- 
tial). This allows us to discuss the experimental observ- 
ability of these corrections to mean-field theory. 

The Hamiltonian of A^ interacting particles in an 
isotropic harmonic trap with frequency uj is given by 
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where V is the interatomic (pseudo-)poteiitial between 
two particles. A hard-core interaction with radius a is 
customary. 

The partition function Z of the system with inverse 
temperature /3 = [ksT)^^ is given by the trace of the 
symmetrized density matrix p = e~^^ over all states. Z 
satisfies the usual convolution equation: 
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Here r = P/M, R is the 3N-dimensional vector 
of the particle coordinates R — (ri, ^2, r^r), and 
R^ denotes the vector with permuted labels: R^ = 
(rp(i-), rp(2) , rp(^)) [|l5|. As explained elsewhere ||] 
the QMC calculation relies on virtually exact formulas 
for the density matrices p{R, R' , t) at the higher temper- 
ature 1/r and performs the integral over i?, i?2 , . . • , Rm 
as well as the sum over all permutations P in Eq. (||) by 
Monte Carlo sampling. Special data-handling techniques 
allow to cope with very large atom numbers N. 

As in |9) we consider a model system of 10,000 parti- 
cles, correponding to a critical temperature of ksT^ — 
20.25?ia; and a hard-core potential with radius a = 
0.0043ao (flo = (h/mujy/^). These values are typical for 
most *^Rb experiments which are now in operation. We 
generally perform computations at different values of r, 
and extrapolate to r ^ 0, in which limit the QMC formu- 
las for the density matrices in Eq. (|^) become manifestly 
exact. 

In BEC, a single quantum state is occupied by a very 
large number of bosons. In the presence of interactions 
this state is a complicated many-body wave function. As 
far as one particle properties are concerned it can, how- 
ever, be described by a single particle wave function ipoi''^) 
determined by the modified Gross-Pitaevskii equation 



2m 



-^r' + U [no(r) + 2nT(r)] - poj M^) 



- 



(3) 



where the densities of condensed no(r) = Nq x ji/'olr)!^ 
and thermally excited particles riT(r) account for inter- 
actions between these particles. The strength of the in- 
teraction is given by U — ATrh^a/m, where a denotes 
again the hard sphere radius. The factor of 2 in front of 
riT accounts for quantum statistical exchange energy. In 
the semiclassical HF approximation the thermal density 
is given by 
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-(muj^r^ /2+2U n{r)-fiT)/kBT 



Here, the thermal wavelength At = fiy/^Ti /mkBT and 
the Bose function g2,/2{z) = Y^=\Z^ I f"^"^ have been 
used. In order to obtain the density distributions of the 
condensate and of the thermal component, Eqs (|^) and 
(^) have to be solved conjointly at the same chemical 
potential p, = pa = pT with the constraint of the fixed 
particle number N (different notations for the chemical 
potential appear for later convenience, cf. below). HF 
neglects the collective excitations of a more fundamen- 
tal Bogoliubov theory. However, it was observed that 
these excitations do not contribute significantly to ther- 
modynamic properties at temperatures well above the 
chemical potential jl^. On the other hand, collective 
excitations are fully included in our QMC calculation, as 
recently confirmed by a comparison of two-particle cor- 
relation functions . 

The precise solution of the HF equation is quite diffi- 
cult. We have obtained identical results both with an it- 
erative procedure and an interpolation-minimization 
routine cf. j^. In a small window around the transition 
temperature, the HF equations have no solution. 
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FIG. 1. Number density N{r) from HF and the exact 
QMC calculations compared to the non-interacting Bose gas, 
both for = 10, 000 particles in an isotropic trap with 
l3 = 0.07/hLu (T~0.7r°). 



The raw output of our Monte Carlo simulations con- 
sists in histograms of N{r), the number of particles at a 
distance r from the center of the trap. In Fig. 1, we com- 
pare the QMC results for N{r) with the HF solution far 
below the critical temperature. We notice an excellent 
agreement between the two approaches. For comparison 
we also show the result for 10, 000 non-interacting bosons 
at the same parameters. The difference is very impor- 
tant. In Fig. 2, we show corresponding data for a tem- 
perature closer to the critical point. Important qualita- 
tive information can be obtained directly from Fig. 1 and 
Fig. 2 without further data analysis i) with respect 
to the non-interacting solution, the number of condensed 
particles (both for HF and QMC) is decreased, while their 
distribution is widened. This can be seen from the first 
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peak in N{r), which is shifted to larger values of r. Notice 
that, even far away from the center of the trap, the dis- 
tribution of the interacting gas is quantitatively different 
from the ideal gas. We have checked that a temperature 
determination of these atoms within the ideal gas model 
fails by about 5%. ii) With respect to the mean- field so- 
lution, the number of condensed particles in the exact nu- 
merical calculation is mcreased. This second point nicely 
checks with recent QMC and renormalization group 
calculations jl3| for the isotropic case. These calcula- 
tions indicate that the degeneracy parameter is changed 
by interaction effects beyond mean-field and that for a 
dilute gas the critical temperature is increased linearly 
in a/ At. It might be suspected that this increase is a 
mere consequence of quasi-particle excitations. However, 
it has been shown within the Popov approximation that 
quasi-particles lead rather to a minute decrease of the 
critical temperature pH . 



N(r) 
0.12 




0.0 5.0 10.0 15.0 

FIG. 2. N{r) from HF and the exact QMC calculations 
for the same system as in Fig. 1, but at a higher temperature, 
13 = 0.056/;itJ (T ~ 0.88 T°). The inset shows the two curves 
in the overlap region. 

To analyze the data we have extracted the number of 
condensed particles as well as their distribution function 
from the QMC calculation. In Ref. Q both were ob- 
tained from the permutation-cycle lengths (cf . also ) . 
Here, however, we have preferred to directly "fit" the 
QMC-generated histograms of N{r) from a very large 
number of samples to smooth curves rio(r) and nT{r), 
as has been done to analyze experimental in situ mea- 
surements. We have found a nearly perfect fit to the 
data, on the local density approximation level, by using 
functional forms for no(r) and nxir) as given in Eqs (j^) 
and with unrestricted 'chemical potentials' /io and 
fj,T- Eqs (^) and (^) are used, but the fit procedure is 
not a HF approximation in disguise. This can be seen 
by a simple counting argument: In the HF calculation, 
the single parameter fj, allows to satisfy the constraint on 
the total number of particles. The HF solution, when- 
ever it exists, usually is unique. In the present case. 



under the same conditions, we have one more parameter, 
which we fix by the condition of minimizing the mean- 
square displacement between the data points and the in- 
terpolating function. The obtained was compatible 
with purely statistical deviations for histograms contain- 
ing more than 10^ data points. We note in passing that 
the density p(r) = N{r)/AiTr'^ and the component densi- 
ties arc obtained with very good precision after perform- 
ing the fit of the raw-data histogram rather than from a 
direct rescaling of the data. 

This fit could also provide a scheme for the accurate 
determination of T within experimental in-situ measure- 
ments. In current experiments the temperature is mostly 
determined by a time-of-flight method. Due to collisions 
during the expansion and the non-instantaneous switch- 
off of the trap this method limits the temperature deter- 
mination to 5%, so that a more accurate determination 
of the critical temperature and the condensate fraction 
has not yet been possible pO| , pl| . 

The use of the fit allows us to compute N() vs. tempera- 
ture for the set of parameters chosen. The results, extrap- 
olated to T — *■ 0, are plotted in Fig. 3 together with the 
HF result, which gives a consistently smaller value of Ni^. 
It should be noticed that the HF equations, as mentioned, 
have no solution for a range of temperatures close to Te- 
la fact, on approaching this window from below the effec- 
tive chemical potential /Zgg = —mu)'^r'^/2 — 2 Un{r) + fiQ 
in Eq. approaches zero (cf. inset of Fig. 3) in the 
overlap region of the two components. Since the Bose 
function (73/2(2) diverges for z > 1, this causes a discon- 
tinuity of the HF solution. This problem is a finite size 
effect which is caused by the nonvanishing kinetic energy 
of the condensate wave function. 
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FIG. 3. Extrapolated condensate fraction No/N from HF 
(lower) and from QMC (upper) vs T in units of hui. The grey 
line interpolates the QMC data. The inset illustrates the in- 
stability of the HF solution as the effective chemical potential 
|^gg-(r)| approaches zero. 

The correction to mean-field theory shown in Fig. 3 
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provides definite support for the scenario found in the 
spatially homogeneous Bose gas: Interactions lead to an 
increased tendency towards Bose condensation ]l2| , p^ . 
There is however an important difference of context: in 
a homogeneous system the density at each point is un- 
changed by interactions and the shift in Tc corresponds 
directly to corrections to HF. As we have seen in the 
present case the trap is much more complicated since the 
finite-size effects, the mean-field itself, and finally the 
corrections to mean-field all influence the density pro- 
file and the condensate fraction. However, the finite- 
size effects and the mean-field itself are already included 
in HF. Therefore the difference between QMC and HF 
become directly comparable to the homogeneous sys- 
tem. In Ref. the following behavior is proposed: 
6Tc/T° = (Tc - T°)/T° ~ 0.34 ap°-^^ ~ 0.46a/ATO. In 
our case we would expect 6Tc ~ 0.065. This corresponds 
roughly to the apparently constant offset in T of the two 
lines in Fig. 3. 

It has been suggested @ that a direct measurement 
of the density in the center of the trap p{r — 0) would be 
useful to determine the deviation of the degeneracy pa- 
rameter from the mean-field result, p{0)X^ = 2.61. Even 
though we are able to extract p{0) from QMC or from ex- 
perimental data very precisely, the onset of condensation 
is not very sharp and the finite-size effects prevent a pre- 
cise determination of the critical degeneracy parameter 
p2[ . We rather advocate a direct comparison between 
HF (which is available for finite systems) and the data. 

It is evident from an experimental point of view that 
the observation of these corrections necessitates a high 
resolution for T and N. Even for ideal experiments the 
intrinsic sample-to-sample fluctuations may prevent their 
observability. Within QMC, we have determined these 
fluctuations from repeated measurements of the complete 
distribution of one whole sample. To be able to dis- 
tinguish clearly the mean particle distribution from the 
mean-field result, one single configuration of the 10, 000 
particles is not sufficient. We find that an average of at 
least 50 independent samples has to be taken in order to 
detect that there is a difference between HF and the ex- 
act result in the sense of a Kolmogorov-Smirnov test [p3| . 
For a precise evaluation (as in Fig. 3) a much larger data 
pool is necessary. Experimentally this seems to require a 
nondestructive measurement scheme. 

In conclusion, we are convinced that comparison be- 
tween theory and experiment on the level treated in this 
paper is conceivable. An experimental access to the 
mean-field corrections will be of outstanding theoretical 
interest [|4|. 
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